﻿module dnaTools

open System

// 4.1 - Meets C1 - C6
let convertAT n l k1 k4 = 
    let M = DWFx.DetWords n l k1 k4 DWFx.BIN_GLUE_SET
    for i = 0 to M.GetLength(0) - 1 do
        for j = 0 to M.GetLength(1) - 1 do
            if (M.[i, j] = '0') then
                M.[i, j] <- 'A'
            else
                M.[i, j] <- 'T'
        done
    done
    M


let outputAT n l k1 k4 = 
    let M = convertAT n l k1 k4
    let bob = M.GetLength(1)
    let Mbox = Array.init bob (fun i -> utils.getColumn i M)
    tester.testBasicHammingConstraint Mbox k1
    tester.testShiftingHammingConstraint Mbox k4
    Mbox


let appendC (M : char [,]) k2 k3 k5 k6 = 
    let mutable k = 0
    for i = 0 to 3 do
        if (max k2 k3) >= (max k5 k6) then
            k <- max k2 k3
        elif (max k3 k5) >= (max k2 k6) then
            k <- max k3 k5
        elif (max k5 k6) >= (max k2 k3) then
            k <- max k5 k6
    done

    // n is supposed to be 1, l is supposed to be 0

    let n = M.GetLength(1)
    let l = M.GetLength(0)
    let mutable combinedArray = Array2D.init (l+k) n (fun i j -> 'B')

    // manually append the C character for the n by k of the matrix
    // fill in only DetWords matrix for the rest
    for p = 0 to (n-1) do
        let mutable j = 0
        while j < (k-1) do
            for j = 0 to (k-1) do
                combinedArray.[j,p] <- 'C'
            done
            j <- j+1
        done
        while j < (l+k-1) do
            for j = k to (l+k-1) do
                combinedArray.[j,p] <- M.[(j-k),p]
            done
            j <- j+1
        done
    done
    combinedArray

// Satisfies contraint C1, C2, C3, C4, C5, and C6
let FastDWD1to6 n el k1 k2 k3 k4 k5 k6 =
    appendC (convertAT n el k1 k4) k2 k3 k5 k6


let outputCA n l k1 k2 k3 k4 k5 k6 =
    //let M = convertAT n l k1 k4
    //let final_M = appendC M k2 k3 k5 k6
    let final_M = FastDWD1to6 n l k1 k2 k3 k4 k5 k6
    let bob = final_M.GetLength(1)
    let Mbox = Array.init bob (fun i -> utils.getColumn i final_M)
    tester.testBasicHammingConstraint Mbox k1
    tester.testReverseComplementaryConstraint Mbox k2
    tester.testSelfReverseComplementaryConstraint Mbox k3
    tester.testShiftingHammingConstraint Mbox k4
    tester.testShiftingReverseComplementaryConstraint Mbox k5
    tester.testShiftingSelfReverseComplementaryConstraint Mbox k6
    Mbox

    (*
    let createTable =
        let p = final_M.GetLength(0)
        let q = final_M.GetLength(1)
        let mutable tableString = ""
        for i = 0 to (p-1) do
            tableString <- tableString + "<tr><td>Word " + string (i+1) + ": &nbsp; 3'</td>"
            for j = 0 to (q-1) do
                if final_M.[i,j] = 'A' then
                    tableString <- tableString + "<td style=\"background-color: #5050ff; padding: .3em; color: black;\">" + "A" + "</td>"
                elif final_M.[i,j] = 'C' then
                    tableString <- tableString + "<td style=\"background-color: #e00000; padding: .3em; color: black;\">" + "C" + "</td>"
                elif final_M.[i,j] = 'T' then
                    tableString <- tableString + "<td style=\"background-color: #e6e600; padding: .3em; color: black;\">" + "T" + "</td>"
                elif final_M.[i,j] = 'G' then
                    tableString <- tableString + "<td style=\"background-color: #00c000; padding: .3em; color: black;\">" + "G" + "</td>"
                //tableString <- tableString + "<td>" + final_M.[i,j].ToString() + "</td>"
            done
            tableString <- tableString + "<td>5'</td></tr><tr><td>&nbsp;</td></tr>"
        done
        tableString
            
    sprintf "%s" createTable
    *)



// 4.2 - Meets C1 through C7
let appendK n l k1 k2 k3 k4 k5 k6 = 
    let mutable k = 0
    for i = 0 to 3 do
        if (max k2 k3) >= (max k5 k6) then
            k <- max k2 k3
        elif (max k3 k5) >= (max k2 k6) then
            k <- max k3 k5
        elif (max k5 k6) >= (max k2 k3) then
            k <- max k5 k6
    done
    let B = DWFx.DetWords n l k1 k4 DWFx.BIN_GLUE_SET
    let new_B = Array2D.init (l+k+k) n (fun i j -> '0')
    for p = 0 to (n-1) do
        let mutable j = 0
        while j < (k-1) do
            for j = 0 to (k-1) do
                new_B.[j,p] <- '1'
            done
            j <- j+1
        done
        while j < (l+k-1) do
            for j = k to (l+k-1) do
                new_B.[j,p] <- B.[(j-k),p]
            done
            j <- j+1
        done
        while j < (l+k+k-1) do
            for j = l+k to (l+k+k-1) do
                new_B.[j,p] <- '1'
            done
            j <- j+1
        done
    done
    new_B

let convertACG (B : char [,]) (y : float) =
    let n = B.GetLength(1)
    let el = B.GetLength(0)
    let lengthY = int (ceil (float el*y))
    // let modBY = l % y
    let listY = Array.init lengthY (fun i -> int(y*(float(i+1)))) // GENERATE MY LIST!!!
    //for i = 0 to listY.Length-1 do
    //    listY.[i] <- int y*(i+1)
    for j = 0 to n-1 do

        // fill in the y positions // I TOOK THE i-1 out to make i
        for i in listY do
            if B.[i,j] = '0' then
                B.[i,j] <- 'C'
            elif B.[i,j] = '1' then
                B.[i,j] <- 'G'
        done   

        // fill in the rest
        for p = 0 to el-1 do
            if B.[p,j] = '0' then
                B.[p,j] <- 'A'
            elif B.[p,j] = '1' then
                B.[p,j] <- 'T'
        done
    done
    B

let outputB n l k1 k2 k3 k4 k5 k6 y =
    let B = appendK n l k1 k2 k3 k4 k5 k6
    let final_B = convertACG B y
    let bob = final_B.GetLength(1)
    let Mbox = Array.init bob (fun i -> utils.getColumn i final_B)
    tester.testBasicHammingConstraint Mbox k1
    tester.testReverseComplementaryConstraint Mbox k2
    tester.testSelfReverseComplementaryConstraint Mbox k3
    tester.testShiftingHammingConstraint Mbox k4
    tester.testShiftingReverseComplementaryConstraint Mbox k5
    tester.testShiftingSelfReverseComplementaryConstraint Mbox k6
    Mbox
    (*
    let createTable =
        let p = final_B.GetLength(0)
        let q = final_B.GetLength(1)
        let mutable tableString = ""
        for i = 0 to (p-1) do
            tableString <- tableString + "<tr><td>Word " + string (i+1) + ": &nbsp; 3'</td>"
            for j = 0 to (q-1) do
                if final_B.[i,j] = 'A' then
                    tableString <- tableString + "<td style=\"background-color: #5050ff; padding: .3em; color: black;\">" + "A" + "</td>"
                elif final_B.[i,j] = 'C' then
                    tableString <- tableString + "<td style=\"background-color: #e00000; padding: .3em; color: black;\">" + "C" + "</td>"
                elif final_B.[i,j] = 'T' then
                    tableString <- tableString + "<td style=\"background-color: #e6e600; padding: .3em; color: black;\">" + "T" + "</td>"
                elif final_B.[i,j] = 'G' then
                    tableString <- tableString + "<td style=\"background-color: #00c000; padding: .3em; color: black;\">" + "G" + "</td>"
                //tableString <- tableString + "<td>" + final_M.[i,j].ToString() + "</td>"
            done
            tableString <- tableString + "<td>5'</td></tr><tr><td>&nbsp;</td></tr>"
        done
        tableString
            
    sprintf "%s" createTable
    *)

let outputBH n l k1 k2 k3 k4 k5 k6 y =
    let B = appendK n l k1 k2 k3 k4 k5 k6
    let final_B = convertACG B y
    final_B

// 4.4 Meets C1 through C8 and uses BreakRuns
let transformC8A (B : char [,]) k2 k3 k4 k5 k6 y (d : int) =
    let n = B.GetLength(0) // length of glues
    let l = B.GetLength(1) // number of glues
    // SECOND coordinate is the x (it's flipped
    let mutable k = 0
    // perhaps add k4 in here, but k4 is already factored into DetWords
    for i = 0 to 3 do
        if (max k2 k3) >= (max k5 k6) then
            k <- max k2 k3
        elif (max k3 k5) >= (max k2 k6) then
            k <- max k3 k5
        elif (max k5 k6) >= (max k2 k3) then
            k <- max k5 k6
    done
    let mutable appB = Array2D.init (n+k+k) l (fun i j -> '0')
    let append1 =
        for i = 0 to (l-1) do // num of glues
            for j = 0 to (n-1) do // just a small length
                appB.[j,i] <- '1'
            done
            for j = k to (n+k-1) do // the appropriate length
                appB.[j,i] <- B.[j-k,i]
            done
            for j = k+l to (n+k+k-1) do
                appB.[j,i] <- '1'
            done
        done
        appB

    // Step 2 of algorithm
    // set up variables
    let el = appB.GetLength(0)
    let new_l = int (floor (float (el/d)))
    let listDiv = Array.init new_l (fun i -> 0)
    // set array sizes here (Z-0, Z-1, etc)
    for i = 0 to listDiv.Length-1 do
        listDiv.[i] <- int d*(i+1)
    done
    // start populating the matrices, we go down each column 
    // until the "end" of the matrix, rather than across the row

    let mutable p = 0
    // if length directly divides into d
    if el % d = 0 then
        for j = 0 to listDiv.Length-1 do
            // even indices
            if listDiv.[j] % 2 = 0 then
                for q = p to listDiv.[j] - 1 do
                    for i = 0 to (l-1) do
                        if appB.[q,i] = '0' then
                            appB.[q,i] <- 'A'
                        elif appB.[q,i] = '1' then
                            appB.[q,i] <- 'T'
                    done
                    p <- p + 1
                done
            elif listDiv.[j] % 2 = 1 then
                for q = p to listDiv.[j] - 1 do
                    for i = 0 to (l-1) do
                        if appB.[q,i] = '0' then
                            appB.[q,i] <- 'C'
                        elif appB.[q,i] = '1' then
                            appB.[q,i] <- 'G'
                    done
                    p <- p + 1
                done
        done

    // if length does not directly divide into d
    else
        // remainder will be our last "matrix set"
        let rem = int (el % d)
        let newlistDiv = Array.init (new_l+1) (fun i -> 0)
        // set array sizes here (Z-0, Z-1, etc)
        for i = 0 to newlistDiv.Length-2 do
            newlistDiv.[i] <- int d*(i+1)
        done
        newlistDiv.[newlistDiv.Length-1] <- rem
        for j = 0 to listDiv.Length-1 do
            // even indices
            if listDiv.[j] % 2 = 0 then
                for q = p to listDiv.[j] - 1 do
                    for i = 0 to (l-1) do
                        if appB.[q,i] = '0' then
                            appB.[q,i] <- 'A'
                        elif appB.[q,i] = '1' then
                            appB.[q,i] <- 'T'
                    done
                    p <- p + 1
                done
            elif listDiv.[j] % 2 = 1 then
                for q = p to listDiv.[j] - 1 do
                    for i = 0 to (l-1) do
                        if appB.[q,i] = '0' then
                            appB.[q,i] <- 'C'
                        elif appB.[q,i] = '1' then
                            appB.[q,i] <- 'G'
                    done
                    p <- p + 1
                done
        done
    appB

// 4.4 part B (Lemma 20) uses BRgamma, not normal gamma (thus is a float, not a int)
//    In order to accomplish this lemma, we must make four different matrices for each sub-word.
//    Thus, we first create B0, where a subword is of length (d-1) except for the remainder.
//    For each subword, the opposite bit is appended to the right. This will create B1.
//    For B2, each word is appended by 1 on the left and 0 on the right.
//    For B3, create an array that is length-d full of 1s ending in 0. Add to left and right.
//    For the final one, for the characters up to ceil(gamma*l) change 0 to C, 1 to G.
//    Everything else should be 0 to A and 1 to T
let transformC8B (B : char [,]) k2 k3 k4 k5 k6 (gamma : float) d =
    let n = B.GetLength(1)
    let mutable k = 0
    for i = 0 to 3 do
        if (max k2 k3) >= (max k5 k6) then
            k <- max k2 k3
        elif (max k3 k5) >= (max k2 k6) then
            k <- max k3 k5
        elif (max k5 k6) >= (max k2 k3) then
            k <- max k5 k6
    done
    let l0 = B.GetLength(0)
    // Run Step 1
    let numSWords = int(ceil (float ((float l0)/(float (d - 1)))))
    let l1 = numSWords + l0
    let B1 = Array2D.init l1 n (fun i j -> '0')
    let mutable listH = List.init (numSWords+1) (fun i -> if i = 0 then 0 else (d-1)*(i)-1) // index where each subword ends
    for i = 0 to n-1 do // for each word i
        let mutable h = 0 // some counter h
        for q = 0 to numSWords-1 do // for each character q of a word up to some limit
            // copy the old values in first (d-1) is the subword length multiplied by which subword it is
            for j = listH.[q] to listH.[q+1] do // for each character j from h to (d-1)*(q+1) - 1
                // put the orig B values here
                B1.[j,i] <- B.[j-(q),i]
            done
            if B1.[listH.[q+1],i] = '1' then
                B1.[listH.[q+1]+1,i] <- '0'
            else
                B1.[listH.[q+1]+1,i] <- '1'
            h <- h + 1
        done
    done

    // Run Step 2
    let l2 = l1 + 2
    let B2 = Array2D.init l2 n (fun i j -> '0')
    for i = 0 to n-1 do
        B2.[0,i] <- '1'
        for j = 1 to l2-2 do
            B2.[j,i] <- B1.[j-1,i]
        done
        B2.[l2-1,i] <- '0'
    done

    // Run Step 3
    let numS2 = int(ceil (float (k/(d-2))))
    let l3 = l2 + 2*numS2*d
    let B3 = Array2D.init l3 n (fun i j -> '0')
    let appendMe = Array.init d (fun i -> '1')
    appendMe.[d-1] <- '0'
    for i = 0 to n-1 do
        // fill in the left appendix
        for j = 1 to numS2 do
            for p = 0 to d-1 do
                B3.[p*j,i] <- appendMe.[p]
            done
        done
        // fill in the middle
        for j = numS2*d to l3-numS2*d-1 do
            B3.[j,i] <- B2.[j-numS2*d,i]
        done
        // fill in the right appendix
        for j = 1 to numS2 do
            for p = 0 to d-1 do
                B3.[p*j+l2+numS2*d,i] <- appendMe.[p]
            done
        done
    done

    // Run Step 4
    let gammaLength = int (ceil (gamma*(float l3)))
    for i = 0 to n-1 do
        for j = 0 to gammaLength-1 do
            if B3.[j,i] = '0' then
                B3.[j,i] <- 'C'
            else
                B3.[j,i] <- 'G'
        done
        for j = gammaLength to l3-1 do
            if B3.[j,i] = '0' then
                B3.[j,i] <- 'A'
            else
                B3.[j,i] <- 'T'
        done
    B3

let outputC8BArray n l k1 k2 k3 k4 k5 k6 (gamma:float) d =
    let B = DWFx.DetWords n l k1 k4 DWFx.BIN_GLUE_SET
    let finalB = transformC8B B k2 k3 k4 k5 k6 gamma d
    let bob = finalB.GetLength(1)
    let Mbox = Array.init bob (fun i -> utils.getColumn i finalB)
    printfn "M: %A" finalB
    printfn "mbox: %A" Mbox
    tester.testBasicHammingConstraint Mbox k1
    tester.testReverseComplementaryConstraint Mbox k2
    tester.testSelfReverseComplementaryConstraint Mbox k3
    tester.testShiftingHammingConstraint Mbox k4
    tester.testShiftingReverseComplementaryConstraint Mbox k5
    tester.testShiftingSelfReverseComplementaryConstraint Mbox k6
    tester.testConsecutiveBaseConstraint Mbox d
    tester.testGCContentConstraint Mbox gamma
    finalB

let outputC8B n l k1 k2 k3 k4 k5 k6 (gamma:float) d =
    let B = DWFx.DetWords n l k1 k4 DWFx.BIN_GLUE_SET
    let finalB = transformC8B B k2 k3 k4 k5 k6 gamma d
    let bob = finalB.GetLength(1)
    let Mbox = Array.init bob (fun i -> utils.getColumn i finalB)
    tester.testBasicHammingConstraint Mbox k1
    tester.testReverseComplementaryConstraint Mbox k2
    tester.testSelfReverseComplementaryConstraint Mbox k3
    tester.testShiftingHammingConstraint Mbox k4
    tester.testShiftingReverseComplementaryConstraint Mbox k5
    tester.testShiftingSelfReverseComplementaryConstraint Mbox k6
    tester.testGCContentConstraint Mbox gamma
    tester.testConsecutiveBaseConstraint Mbox d
    Mbox
    (*
    let createTable =
        let p = finalB.GetLength(0)
        let q = finalB.GetLength(1)
        let mutable tableString = ""
        for i = 0 to (p-1) do
            tableString <- tableString + "<tr><td>Word " + string (i+1) + ": &nbsp; 3'</td>"
            for j = 0 to (q-1) do
                if finalB.[i,j] = 'A' then
                    tableString <- tableString + "<td style=\"background-color: #5050ff; padding: .3em; color: black;\">" + "A" + "</td>"
                elif finalB.[i,j] = 'C' then
                    tableString <- tableString + "<td style=\"background-color: #e00000; padding: .3em; color: black;\">" + "C" + "</td>"
                elif finalB.[i,j] = 'T' then
                    tableString <- tableString + "<td style=\"background-color: #e6e600; padding: .3em; color: black;\">" + "T" + "</td>"
                elif finalB.[i,j] = 'G' then
                    tableString <- tableString + "<td style=\"background-color: #00c000; padding: .3em; color: black;\">" + "G" + "</td>"
            done
            tableString <- tableString + "<td>5'</td></tr><tr><td>&nbsp;</td></tr>"
        done
        tableString
            
    sprintf "%s" createTable
    *)

// *** TODO: Find a way to display the complementary value so that it will display correctly
// *** MAYBE WE JUST PASS ON THIS ARRAY FOR FURTHER PROCESSING?
// outputs a usable array version of before
let outputC8AArray n l k1 k2 k3 k4 k5 k6 y d =
    if l<=0 then Array2D.zeroCreate 0 0
    else
        printf "\tGenerating DetWords Matrix...\n"
        let B = DWFx.DetWords n l k1 k4 DWFx.BIN_GLUE_SET
        //printf "result: %A\n" B
        printf "\tGenerating Output Array...\n"
        let finalB = transformC8A B k2 k3 k4 k5 k6 y d
        let bob = finalB.GetLength(1)
        let Mbox = Array.init bob (fun i -> utils.getColumn i finalB)
        printfn "M: %A" finalB
        printfn "mbox: %A" Mbox
        tester.testBasicHammingConstraint Mbox k1
        tester.testReverseComplementaryConstraint Mbox k2
        tester.testSelfReverseComplementaryConstraint Mbox k3
        tester.testShiftingHammingConstraint Mbox k4
        tester.testShiftingReverseComplementaryConstraint Mbox k5
        tester.testShiftingSelfReverseComplementaryConstraint Mbox k6
        tester.testConsecutiveBaseConstraint Mbox d
        finalB

let outputC8A n l k1 k2 k3 k4 k5 k6 (y:float) d =
    let B = DWFx.DetWords n l k1 k4 DWFx.BIN_GLUE_SET
    let finalB = transformC8A B k2 k3 k4 k5 k6 y d
    let bob = finalB.GetLength(1)
    let Mbox = Array.init bob (fun i -> utils.getColumn i finalB)
    tester.testBasicHammingConstraint Mbox k1
    tester.testReverseComplementaryConstraint Mbox k2
    tester.testSelfReverseComplementaryConstraint Mbox k3
    tester.testShiftingHammingConstraint Mbox k4
    tester.testShiftingReverseComplementaryConstraint Mbox k5
    tester.testShiftingSelfReverseComplementaryConstraint Mbox k6
    tester.testConsecutiveBaseConstraint Mbox d
    Mbox